Richardson's pair diffusion and the stagnation point structure of turbulence 
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DNS and laboratory experiments show that the spatial distribution of straining stagnation points 
in homogeneous isotropic 3D turbulence has a fractal structure with dimension D s = 2. In Kinematic 
Simulations the time exponent 7 in Richardson's law and the fractal dimension D s are related by 
7 = G/D s . The Richardson constant is found to be an increasing function of the number of straining 
stagnation points in agreement with pair duffusion occuring in bursts when pairs meet such points 
in the flow. 
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The rate with which pairs of points separate in phase 
space or in physical space is of central importance to 
the study of dynamical systems. Pairs of points in the 
phase space of a low-dimensional chaotic dynamical sys- 
tem separate exponentially. In chaotic advection pairs 
also separate exponentially |lj leading to exponentially 
fast stirring and a high potential for mixing. In turbulent 
flows, however, pairs of fluid elements separate on average 
algebraically [g, Q . Turbulent flows have a very wide 
range of excited length- and time-scales and are therefore 
fundamentally different from both low-dimensional dy- 
namical systems and chaotic advection flows. Attempts 
have been made to model the relative diffusion of fluid 
element pairs in terms of Langevin type equations based 
on the assumption that relative Lagrangian accelerations 
are Markovian in time [|. These models can repro- 
duce the right algebraic time growth of separation statis- 
tics of fluid elements in turbulent flows; but they fail to 
explain the very large values taken by the flatness fac- 
tor of Lagrangian relative velocities in Direct Numerical 
Simulations of isotropic turbulence || 0, § ■ In fact these 
models based on Markovian acceleration statistics under- 
estimate this flatness factor by as much as one order of 
magnitude. 

The algebraic growth of Lagrangian separation statis- 
tics can also be accurately reproduced by Kinematic Sim- 
ulations §,@0. These are models of turbulent diffu- 
sion based on kinematically simulated turbulent velocity 
fields which are non-Markovian (not delta-correlated in 
time), incompressible and consistent with up to second 
order statistics of the turbulence such as energy spec- 
tra. Kinematic Simulations are interesting in particular 
because they do reproduce the very high flatness fac- 
tors of Lagrangian relative velocities [||. The mecha- 
nism by which fluid element pairs separate in Kinematic 
Simulations (KS) might therefore be comparable to the 
one in turbulent flows and is clearly different from the 
Wiener process which causes fluid element pairs to sepa- 
rate in Lagrangian models of relative diffusion based on 
Langevin type equations. But what is this mechanism 
and why can it give rise to the algebraic growth of rela- 



tive separations? 

Richardson's law of turbulent relative diffusion states 
that the mean square distance between two fluid elements 
A 2 is proportional to the third power of time t, i.e. 



A 2 (t) 



(1) 



where 7 = 3, G is a universal dimcnsionlcss constant and 
L and v! are the integral length-scale and the rms ve- 
locity of the turbulence respectively. Richardson's law 
is expected to be valid in homogeneous isotropic turbu- 
lence and for times t such th at A 2 is within the inertial 
range of scales, i.e. i) 2 C A ! < L 2 , where 77 is the Kol- 
mogorov microscale. Fluid element pairs follow close tra- 
jectories for long stretches of time and separate violently 
when they meet straining flow regions around stagnation 
points (straining stagnation points) || O, [fj, |l3). || 
noted that the very high flatness factors of Lagrangian 
relative velocities are consistent with this conjecture. 

[]ll"| found evidence of a fractal spatial distribution of 
straining flow regions in KS homogeneous isotropic tur- 
bulent velocity fields. In this paper we quantify their ob- 
servation by showing that the number of straining stag- 
nation points per unit volume in KS, Direct Numerical 
Simulations (DNS) and laboratory experiments of homo- 
geneous isotropic turbulence is given by 



C S L- 



D„ 



(2) 



where D s = 2 and C s is a dimensionless number. The 
exponent D s can be interpreted as a fractal dimension. 
Both (0) and (|) hold for L/77 > 1. 

We now describe the DNS and the grid turbulence mea- 
surements used to establish (||) with D s =2. We have used 
DNS data of non-decaying homogeneous isotropic incom- 
pressible turbulence generated by a standard pseudo- 
spectral code with grid resolution of about 2ij and have 
computed the number of stagnation points in instanta- 
neous velocity fields u = u(x) = (u(x), v(x), u>(x)) for a 
variety of Taylor microscale Reynolds numbers Re\ rang- 
ing from 34 to 130. In this paper we focus our interest on 
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FIG. 1: Imaginary versus real part of the complex eigenvalues 
of the gradient of velocity matrix with positive imaginary part 
for 26726 stagnation points in a KS field with L/r/ « 45. Very 
few complex eigenvalues with small (absolute) real part can be 
found. Similar results (with poorer resolutions) were obtained 
with DNS for similar values of L/r]. 



the straining flow regions around stagnation points of u. 
A stagnation point does indeed correspond to a straining 
flow region when the eigenvalues of the velocity gradient 
matrix at this point are all non-zero (see Ottino 1989); 
these are regions where the flow is always straining and 
may or may not be rotating as well. 

We use the Newton-Raphson method (tested against 
the amoeba method in various planar flows) to find all 
stagnation points [|l4|. This is an iterative method and 
requires starting points, which have been taken all over 
the DNS field at points separated by a distance smaller 
than r\. Irrespective of Reynolds number, the vast major- 
ity of stagnation points have turned out to be in straining 
regions. Figure [I] is a scatter plot of the complex eigen- 
values of the gradient of velocity matrix. Notice that 
the probability to find eigenvalues with imaginary part 
much larger than the real part is very low. The only case 
where stagnation points can be non-straining is when two 
eigenvalues are pure imaginary and the third eigenvalue 
is zero. Hence the number of stagnation points is there- 
fore effectively the same as that of straining stagnation 
points, and this number per unit volume is n s . For ev- 
ery Reynolds number, we calculate n s , L and rj and we 
plot n s as a function of L/r) (see figure ||). The relation 
between these two quantities appears to be well fitted by 
(|) with D s = 2. 

Experimental support for (^) with D s =2 has been ob- 
tained from grid generated turbulence in the laboratory. 
The turbulence is homogeneous and isotropic far enough 
from the grid and from the wind tunnel walls. Measure- 
ments of the streamwise air velocity u were taken with a 
hot wire at the centre of the working section of a wind 



° * ' ' 



LA] 



FIG. 2: Number of stagnation points per unit volume versus 
L/r) in DNS (•) and KS for p — 5/3 and different values of V 
(compared to u' = 1): o V= 0, + V= 0.7, v V= 1.0. Dashed 
lines representing n s tx (L/rj) 2 are shown for comparison. 



tunnel at a distance of 50 mesh sizes behind the grid. We 
collected thirty runs of data for thirty different values of 
Re\ ranging from 68 to 130. The sampling frequency was 
30 kHz, enough to resolve the dissipation range except at 
the highest values of Re\ but always enough, however, to 
resolve the Taylor microscale in all our runs. Each data 
set contains more than 100 integral scales. The turbu- 
lence intensity was always smaller than 5% thus allow- 
ing the use of the Taylor hypothesis for the conversion 
of temporal data into spatial data. It is possible, from 
these wind tunnel data, to calculate the number of zero- 
crossings of u per unit length which leaves us with the 
problem of relating this number to the number of stagna- 
tion points per unit volume in homogeneous and isotropic 
turbulence. Each component of u has an instantaneous 
zero-crossing surface in the three-dimensional space of 
the flow. These three instantaneous surfaces may have 
a fractal dimension larger than or equal to 2. Because 
of isotropy, these three fractal dimensions must be the 
same and we denote them by D. Stagnation points of u 
are intersections of the three zero-crossing surfaces. The 
rule of the thumb for estimating the fractal dimension of 
intersections of surfaces is that the codimension is equal 
to the sum of the codimensions of the intersecting sur- 
faces [[l5| . The codimension of each zero-crossing surface 
is D — 2 because surfaces are two-dimensional and the 
codimension of the set of their intersections is D s — be- 
cause points are zero-dimensional. Hence, D s — 3(D — 2). 
We are assuming our DNS observation that the vast ma- 
jority of stagnation points lie in straining regions to be 
also true in grid generated turbulence which is why we 
use the notation D s for the fractal dimension of the stag- 
nation points of u. 

By virtue of the Taylor hypothesis, the zero-crossings 
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of our one-dimensional u data can be viewed as a set 
of point intersections through the zero-crossing surface 
u(x) = 0. From our one-dimensional data the codimen- 
sion D — 2 can be measured from one data set with a 
specific Reynolds number by applying a box-counting al- 
gorithm to the zero-crossings. This method was used by 
Sreenivasan and his colleagues (see |l6| and references 
therein) who found that the fractal codimension D — 2 
is indeed well-defined and equal to 2/3. Another way to 
measure D — 2 is to count the number of zero-crossings 
in different data sets corresponding to different Reynolds 
numbers and relate the zero-crossing numbers per unit 
length to L/t]. However, because we do not resolve r\ in 
some of our runs, we apply a low-pass filter of wavenum- 
ber 2ir/l c in order to remove the poorly resolved dissipa- 
tion range and the high-frequency electronic noise of our 
measuring device. The filter wavenumber is in fact equal 
to 2ir/ri for the smallest Reynolds numbers and turns out 
to be closer to the Taylor microscale wavenumber in most 
cases. This is the method we have applied here and we 
have found that the zero-crossing number per unit length 
is proportional to (L/l c ) 2 / 3 which implies D — 2 — 2/3. 
From D s = 3(D — 2) we deduce D s = 2 in support of our 
DNS findings. 

To explore the relation between the turbulent diffusion 
of fluid clement pairs and the fractal structure of strain- 
ing stagnation points in the flow, i.e. between (|l|) and (||), 
we need to find ways to modify the spatial distribution 
of straining stagnation points and monitor the changes 
in turbulent pair diffusion brought about by such modi- 
fications. Such a study cannot be carried out with cur- 
rent DNS and laboratory experiments of homogeneous 
isotropic incompressible turbulence because the spatial 
distribution of straining stagnation points is determined 
by the Navier-Stokes dynamics and cannot be tampered 
with. KS, however, offers the flexibility to chose the en- 
ergy spectrum at will and thereby modify, as we show 
below, the fractal structure of the set of straining stag- 
nation points. An additional advantage of KS is that the 
Lagrangian pair diffusion statistics it produces compare 
well with DNS results when the energy spectrum chosen 
is that of the DNS turbulence 0. KS also succesfully 
generates [|l2| all the pair diffusion results of the labora- 
tory experiment of |l3] ], 

KS uses turbulent-like velocity fields of the form 

N k 

u = A w A k w cos(k„ • x + uj n t) + 

n=l 

B„ A k„ sin(k„ • x + w n t) (3) 

where Nk (typically of order 100) is the number of modes, 
k„ is a random unit vector (k„ — fc„k„), and the di- 
rections and orientations of A„ and B„ are chosen ran- 
domly under the constraint that they be normal to k ra 
and uncorrelated with the directions and orientations of 
all other wave modes. Note that the velocity field (pi) is 



incompressible by construction, and also statistically sta- 
tionary, homogeneous and isotropic as shown by Q and 
. The amplitudes A n and B n of the vectors A„ and 
B n are determined by A 2 n = B 2 = ^E(k n )Ak n where 
E(k) is the energy spectrum prescribed to be of the form 



E{k) 



u' 2 {p- 1) 
2(L/2tt)p- 1 



(4) 



in the range 2tt/L\ = k\ < k < k]y k = 2ir/rj, and E(k)=0 
otherwise; u' is the rms velocity of the KS turbulent- 
like flow; Ak n = (fc„+i - fe„-i)/2 for 2 < n < N k - 1, 
Afci = (fc 2 - fci)/2 and Ak Nh = (k Nk - fcjv fe -i)/2. The 
distribution of wavenumbers is geometric (see Flohr & 
Vassilicos 2000), specifically k n = fcia™ -1 with a constant 
a determined by Ljr\ — a Nk ~ 1 . The frequencies oj n in (|^) 
are proportional to the eddy-turnover frequency of mode 
n, i.e. uj n = 0.5^k^E(k n ). 

We have varied the power p of the energy spectrum in 
the range 1 < p < 3 without changing u' . For a given 
value of p, we calculate the number of stagnation points 
by the Newton-Raphson method for different values of 
L/rj. It turns out, as in the case of DNS turbulence, 
that the vast majority of stagnation points are straining 
stagnation points irrespective of the values of p and L/rj 
(see figure [l]). We also find that, when p = 5/3, n s is 
proportional to (L/rj) Ds with D s = 2 (see figure ||), in 
agreement with our DNS and wind tunnel results. The 
relation (|J) is found to hold for all values of p between 
1 and 3 in instantaneous realisations of our KS field and 
in fact Ds decreases towards as p increases towards 
3. Varying p in KS is therefore a good way to tamper 
with the fractal structure of the set of straining stagna- 
tion points in the turbulent-like flow and study what the 
effects of this tampering are on turbulent pair-diffusion. 

We have calculated the time-dependence of the mean 
square distance between pairs of fluid elements in KS 
turbulent-like flows for different values of p between 1 
and 2 and have found that (|l|) is valid in the inertial 
range with 



7 = 6/D s 



(5) 



(see figure ||), and that the Richardson constant G is an 
increasing function of D s . For values of p between 2 and 
3 we do not find evidence of a well-defined power law 
(|l|) and of course nothing like (||). The reason, which 
we believe lies behind this pair-diffusion behaviour, is 
that for p between 2 and 3 the power spectrum k 2 E(k) 
of the strain rate field is concentrated at the smaller 
wavenumbers when p > 2 but increases with wavenum- 
ber when p < 2. The energy spectrum of a homogeneous 
and isotropic gaussian velocity field such as for large 
enough Nk is related to the fractal dimension of zero- 
crossing surfaces by Orey's relation p + 2(D — 2) = 3 (see 
p7[) which implies p + 2D s /3 = 3. have shown that, 
as a consequence of pair-diffusion locality, j = 4/(3 — p). 
Hence, ([)]) is effectively a consequence of velocity field 
gaussianity and pair-diffusion locality. 
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FIG. 3: Relation between Richardson's pair diffusion expo- 
nent 7 and the fractal dimension of stagnation points D s . 
The values of 7 are obtained from linear fits of the time 
dependence of A 2 (in log-log scale) over the interval where 
(3??) 2 < A 2 < (L/3) 2 . The length of the error bars is twice 
the r.m.s. of 7 within that interval. In these KS calculations, 
iVfc = 40, L/r] — 10 3 and the initial separation is r//2 for 
2000 different particle pairs. The dashed line shows (pi) for 
comparison. 




FIG. 4: Richardson constant G versus flow average velocity 
V using the parameters of the simulations as in figure B. 



Another way to modify the straining flow structure of 
the KS turbulence is by adding a constant (time and 
space independent) velocity vector V to the KS veloc- 
ity field (g). The KS velocity field defined in (§) is a 
mean-zero velocity field and the addition of V amounts 
to a superposition of a constant mean flow. The aim 
is to look for stagnation points of u + V and monitor 
the changes in the straining stagnation point structure 
caused by changes in V = |V| whilst keeping p constant. 
Our first observation is that ([!]) remains valid with the 



same value of D s but C s decreases with V (see figure 
^). By releasing fluid element pairs in the velocity field 
u + V we can study modifications to the Richardson law 
(1) caused by the mean flow velocity V. The effects on 
(1) parallel those on (g): the exponent 7 remains well- 
defined in the same range of times and is independent of 
V but the Richardson constant G decreases with increas- 
ing V (figure ^). The addition of a constant mean flow 
velocity leaves the scaling exponents 7 and D s of the 
Richardson law and of the fractal straining stagnation 
point structure intact, but reduces the overall number of 
straining stagnation points per unit volume and also the 
overall extent of turbulent pair diffusion. In view of this 
conclusion and also of relation (|J) there is clearly a corre- 
lation between Richardson pair diffusion and the fractal 
spatial distribution of straining stagnation points in the 
turbulcnt-likc flow. Such a correlation hints at a certain 
persistence in time of the streamline structure of the flow 
which causes pairs to separate when they meet straining 
stagnation points. This is consistent with our results that 
G is an increasing function of D s and of 1/V, i.e. an in- 
creasing function of the number of straining stagnation 
points in both cases. 
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